Liquid and crystal phase of dipolar fermions in two dimensions 
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The liquid and crystal phase of a single-component Fermi gas with dipolar interactions are inves- 
tigated using quantum Monte Carlo methods in two spatial dimensions and at zero temperature. 
The dipoles are oriented by an external field perpendicular to the plane of motion, resulting in a 
purely repulsive 1 /r' ! interaction. In the liquid phase we calculate the equation of state as a function 
of the interaction strength and other relevant properties characterizing the Fermi-liquid behavior: 
effective mass, discontinuity at the Fermi surface and pair correlation function. In the high density 
regime we calculate the equation of state of the Wigner crystal phase and the critical density of 
the liquid to solid first order phase transition. Close to the freezing density we also search for the 
existence of a stripe phase, but such a phase is never found to be energetically favorable. 



PACS numbers: 05.30.Fk, 03.75.Hh, 03.75.Ss 

The recent rapid developments in the field of ultracold 
dipolar atoms and molecules has opened up new fasci- 
nating prospects for investigating many-body effects in 
quantum degenerate gases with long-range interactions 
(for a review see, e.g., Ref. pj). In this respect, single- 
layer and multi-layer configurations of two-dimensional 
fermions are particularly intriguing because of the com- 
peting interplay, depending on the strength of the dipolar 
interaction and on the distance between layers, between 
Fermi liquid behavior, supcrfluid pairing, crystal order 
and density- wave instabilities [2h10|. 

Fermionic molecules of 40 K 87 Rb, which can have a 
strong electric dipole moment, have been created using 
coherent transfer of Feshbach molecules to their rovibra- 
tional ground state and have been brought toward 
the quantum degenerate regime Other fermionic 

molecules are now being actively studied experimen- 
tally IH li |- Atomic species with a large magnetic mo- 
ment, such as dysprosium, offer a different possibility of 
realizing degenerate Fermi gases of dipoles that was suc- 
cessfully pursued in the recent experiment reported in 
Ref. Q. 

A particularly simple geometrical arrangement of a 
single-component dipolar Fermi gas in 2D is when the 
dipoles are oriented perpendicular to the plane of motion 
by means of a sufficiently strong external field. This con- 
figuration has been proven to greatly suppress the chem- 
ical reaction rate of molecules, thereby enhancing their 
lifetime flil ]. Here particles at distance r interact via a 
purely repulsive, rotationally symmetric and long range 
1 /r 3 potential. Still the phase diagram at zero tempera- 
ture is expected to be quite rich: interlayer dimers and a 
novel BCS-BEC superfluid crossover are predicted in bi- 
layer systems [f| , while in- plane and out-of-plane density 
ordered phases are predicted in multilayer systems 0, 0] . 
In the case of a single layer, a Fermi liquid with peculiar 
scattering properties is stable at low density [l0| and a 
Wigner crystal emerges at high density, where the clas- 
sical potential energy of dipoles largely exceeds their ki- 
netic energy. For intermediate values of the interaction 



strength an instability at finite wave vector is predicted 
to set in 0, 0, ' driving the system to a stripe phase 
that breaks both rotational and translational symmetry 
(in the direction perpendicular to the stripes) . A similar 
scenario, involving microemulsion phases (e.g. stripes or 
bubbles) is expected for the meting of the Wigner crystal 
at T = in a 2D Coulomb gas[17|. These results are de- 
rived within a mean-field approach: an important ques- 
tion concerns the quantitative determination of the phase 
diagram using more accurate theoretical tools, such as 
quantum Monte Carlo (QMC) techniques [HI ■ 

In this Letter we examine a 2D system of N identi- 
cal fermionic particles of mass m that interact with the 
Hamiltonian 
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where is the distance between particle i and j and d 
is the intensity of the electric (or magnetic) dipole mo- 
ment. The strength of the dipolar interaction is conve- 
niently expressed in terms of the dimensionless parameter 
kprQ, where rg = md 2 /h 2 is the characteristic length of 
the dipole-dipole force and kp = \J Aim is the Fermi wave 
vector of the 2D gas determined by the density n. The en- 
ergy scale set by kp is the Fermi energy tp = Ti 2 k F /2m. 
As a function of kpro we investigate the ground-state 
properties of the Fermi liquid phase including the equa- 
tion of state, the effective mass, the discontinuity of the 
momentum distribution at kp, the pair correlation func- 
tion and the static structure factor. By comparing the 
energy of the Fermi liquid (FL) and of the Wigner crystal 
(WC) phase, we determine the value k F r = 25 ±3 of the 
freezing density. Furthermore, in the region of interac- 
tion strengths close to freezing, we calculate the energy 
corresponding to a stripe phase finding that it is never 
favorable compared to the FL and WC state. The main 
results on the equation of state are shown in Fig. [1] in 
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units of the Hartree-Fock energy 
Ehf 
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corresponding to the lowest order perturbation expansion 
of the FL in the interaction parameter kpr 0, [l^ . 

We use the fixed-node diffusion Monte Carlo method 
(FN-DMC), a projector technique that, starting from an 
antisymmetric trial wave function ip? , finds the state hav- 
ing the lowest energy compatible with the many-body 
nodal surface of ipx that is kept fixed during the calcu- 
lation. The method provides a rigorous upper bound to 
the energy of the fermionic ground-state [19( . 

Simulations are performed in a box of volume V = 
L x L y , where we always take L x < L y . The density 
is n = N/V and we use periodic boundary conditions 
(PBC) in both spatial directions. Since the dipole-dipole 
interaction is long range, the potential energy contribu- 
tion to the Hamiltonian, given by the second term in 
Eq. ([T]), requires a careful treatment 
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where i and j label particles in the simulation cell and the 
vectors r^) + R correspond to the positions of all images 
of particle in the array of replicas of the simulation 
cell. The combination of all images has the same aver- 
age density n of the simulation box and provides a good 
approximation of the homogeneous medium. In the case 
of the Coulomb potential the summation in ([3]) is carried 
out by means of the Ewald method [l9| . For the faster 
convergent 1/r 3 potential the mean interaction energy 
can be evaluated using the simpler formula 



(V) = (V dd ) Rc + v t 



tail 
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where (V dd )R c denotes the sum (J3J) with the constraint 
|rj — r,j — R| < R c and Vtaii = imd 2 / R c is the contribu- 
tion from distances larger than R c assuming a uniform 
distribution of particles. The cut off length R c is chosen 
large enough {R c = 2L X ) to yield an average interaction 
energy (V) that is independent on R c , within statistical 
uncertainty. 

Fermi liquid phase: The trial wave function describing 
the FL phase is assumed of the Jastrow-Slater form 
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Here k Q = (27r/L)«,n») with n%' v = 0,±1,±2,... are 
the wave vectors complying with PBC in the square 
box (L x = L y = L) and f(r) is a non- negative func- 
tion satisfying the boundary condition f'(r = L/2) = 0. 
The short-range behavior of f(r) is of the form f(r) cx 
Kq{2^ r~o/r), where Kq is the modified Bessel function, 



0.9 
0.8 
0.7 
0.6 
0.5 



0.006 
x 0.004 
^ 0.002 
< 
-0.002 



203040 50 6070 
k F r 



0.1 



10 



k F r 



FIG. 1: (color online). Equation of state of the liquid and 
solid phase in units of the Hartree-Fock energy ((2J). Circles 
refer to the liquid and triangles to the solid. The red dashed 
line corresponds to the second-order expansion in Ref. [Io( [. 
The purple dashed horizontal and solid line correspond re- 
spectively to the classical energy of the Wigner crystal and 
to the result of Ref. [28| including the first correction aris- 
ing from the zero-point motion of phonons. Inset: Energy 
difference between the solid and the liquid (blue circles) and 
between the stripe phase and the liquid (black circles). The 
blue solid line is obtained from a best fit to the equation of 
state of the liquid and solid phase. Error bars are smaller 
than the symbols size. 
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FIG. 2: (color online). Finite size effects in a Fermi liquid at 
the density kpro = 0.22. Black symbols and line correspond 
to TABC energies and linear fit. Green symbols and line 
correspond to PBC energies En — (m/?n*)ATjv and linear 
fit. Red symbols correspond instead to PBC energies En — 
ATn- The values extrapolated to the TL are shown with the 
corresponding error bars. 



and fulfills the cusp condition of the dipolc-dipolc poten- 
tial [H. 

A delicate issue related to QMC calculations of the 
equation of state is the extrapolation to the thermody- 
namic limit (TL). In the FL phase apart from the size 
dependence affecting the potential energy contribution, 
which we treat using the procedure in Eq. HI significa- 
tive shell effects are present in the kinetic energy con- 
tribution. We consider closed-shell configurations corre- 
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for which the relative error 
/ jjj tl — 1 1 in the energy of the non- 
interacting gas compared to the TL can be as large as 
~ 1%. An extrapolation method based on FL theory is 
provided by the fitting formula [2lJ 
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which involves the parameter m/m*, determining the in- 
verse effective mass of the particles, and the coefficient 
a of the residual size dependence assumed to be linear 
in 1/iV. Here En is the QMC output energy of the TV- 
particle system with the potential contribution evaluated 
using Eq. [4] The values of En — ATjv are shown as red 
symbols in Fig.O Their scattered dependence on 1/N is 
considerably suppressed if one accounts for the effective 
mass, as it is shown by the green symbols corresponding 
to En — (to/to* )AT/v. A more reliable convergence to 
the TL is obtained using the method of twist-averaged 
boundary conditions (TABC) [Hf. Here the PBC wave 
vectors of the plane waves in the Slater determinant of 
Eq. © are replaced by k Q (0) = {2-k / L){n x a +6 X , K + 6 y ), 
where 9 X , 8 y are continuous variables in the interval [0, 1]. 
In the grand canonical implementation of TABC de- 
scribed in Refs. [23| the wave vectors are constrained 
by |k Q (6>)| < kp and different values of the twist can 
correspond to different numbers of particles. The num- 
ber of particles N and the energy E N are obtained from 
averages over all possible twist angles. With our use of 
TABC we still find a residual size effect AT N [Hj]. The 
extrapolation to the TL is performed using Eq. (|6]) and 
statistical agreement in Etl between PBC and TABC is 
obtained for all values of the density (see Fig. [21). 

The FN-DMC results of the ground-state energy 
are reported in Fig. [T] At low density we find 
good agreement with the result E = Ehf + 
(Nep /2)(kpro) 2 log(1.43fc;?rn), which was derived in 
Ref. [10( including corrections to the lowest order expan- 
sion ©. 

The effective mass to* obtained from Eq. ^ is shown 
in Fig. [3] for different values of kpr - At weak coupling 
our results well reproduce the perturbation expansion re- 
ported in Ref. [ill , but for larger couplings the reduction 
of to* is less pronounced than the perturbative predic- 
tion and m* /to approaches the value 0.4 for kp r$ close 
to freezing. We also calculate the momentum distribu- 
tion n(k) and using the fit = Z8(kp — k) + g(k), 
where 6{x) is the step function and g{k) is a continuous 
function of k, we extract the renormalization factor Z. 
Results are reported in Fig. [3] In the inset of Fig. [3] we 
show nfc in the strongly interacting regime for different 
numbers of particles. It is worth noticing that finite size 
effects on this quantit y ar e much less severe here than in 
the 2D Coulomb gas 25j and are similar to the case of 
hard-disks investigated in Ref. [2(| . 

We finally calculate the pair correlation function g(r) 




FIG. 3: (color online). Effective mass and renormalization 
factor in the liquid phase as a function of the interaction 
strength. The line cor resp onds to the perturbation expan- 
sion for m*/m of Ref. [1Q|. Inset: Momentum distribution 
corresponding to kpro = 20 for different system sizes. 




FIG. 4: (color online). Pair correlation function in the liquid 
and in the crystal phase. The pair correlation function of the 
non-interacting gas is also shown. 




FIG. 5: (color online). Static structure factor in the liquid and 
in the crystal phase. In the liquid phase, solid lines correspond 
to the Fourier transform of g(r) while symbols correspond to 
the direct calculation of S(k). The static structure factor of 
the non-interacting gas is also shown. 
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giving the probability of finding two particles at the dis- 
tance r. Results for different values of k F ro in the FL 
phase are shown in Fig. 01 It is interesting to notice 
that by increasing the interaction strength the short- 
range repulsion increases and a shell structure starts 
to appear on approaching the freezing density. The 
Fourier transform of g(r) yields the static structure fac- 
tor 5(k) = 1 + nj dve iiL ' T [g{r) - 1]. This quantity can 
also be calculated directly in the FN-DMC algorithm by 
evaluating the average of the product of density fluctu- 



ation operators NS{k) = (pk/9-k) 



(E 



Results are reported in Fig. [5] for both estimators [27| ■ 
For large values of k F r$ the direct estimator exhibits a 
more pronounced peak compared to the smoother Fourier 
transform at the wave vector corresponding to the lowest 
non-zero reciprocal lattice vector of the triangular lattice. 

Wigner crystal phase: To describe the WC phase we 
make use of the following trial wave function 



^ T {v u ...,v N ) =n /fai) det[ e -( r - R -) 2 / Q2 ] 



(7) 



density and dipolar strength. This can be understood if 
one considers that the equation of state of the crystal is 
practically independent of statistics while the energy of 
the fermionic fluid is significantly larger than the bosonic 
one. From the equation of state of the FL and WC phase 
in the vicinity of the freezing density one can also esti- 
mate the width of the region where phase separation oc- 
curs driving the first-order liquid to solid transition. By 
imposing equilibrium of pressure and chemical potential 
in the two phases, the coexistence region turns out to 
be d(k F ro) ~ 0.01, a very small value consistent with a 
similar finding in the bosonic case [28| . The pair correla- 
tion function and the static structure factor deep in the 
crystal phase are shown respectively in Fig. |4] and Fig. [5] 
In particular, S(k) exhibits a large peak at k = 1.90fci? 
corresponding to the lowest non-zero wave vector of the 
reciprocal lattice. 

Stripe phase: A pattern of equally spaced stripes is as- 
sumed in the y-direction corresponding to the trial wave 
function 



where the Jastrow correlation term f(r) is the same as in 
the FL phase and the single-particle orbitals in the deter- 
minant are constructed with Gaussians, whose width a 
is a variational parameter, centered at the lattice points 
R TO of the triangular Bravais lattice. In order to enforce 
PBC, both the number of particles TV and the box sizes 



ip T (ri,--. 7 r N ) = Y[f(r lj ) det[e 



(9) 

where y a denotes the y coordinate of the a-th stripe and 
k ax — 2irn ax / L x are the PBC wave vectors in the x- 
direction. The number of fermions is the same in each 
stripe and once multiplied by the number of stripes deter- 



L x and L y must be multiples of a primitive cell, which mines the overall density in the volume V = L x L y . First, 



can be chosen as a rectangle of side lengths l y = \fZl x 
containing two atoms. Extrapolation to the thermody- 
namic limit is obtained using a linear fit in 1/N over the 
FN-DMC energies. The results for the WC equation of 
state are reported in Fig. [TJ It is worth noticing that the 
antisymmetric constraint imposed in the wave function 
(|7J) for particle exchange is negligible for the value of the 
energy. In fact, statistically compatible results are ob- 
tained with a node less wave function of distinguishable 
boltzmanons in agreement with the findings of Ref. pol ] . 
This behavior is expected at large density, where the en- 
ergy of the WC phase is given by the result [28| 
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obtained by including the contribution from the zero- 
point motion of phonons to the purely classical interac- 
tion energy of a system of dipolcs arranged in a trian- 
gular Bravais lattice. The above expansion, holding for 
large k F ro, is shown in Fig. Q] and is indeed approached 
by our QMC results. The difference between the ground 
state energy of the WC and FL phase is shown in the 
inset of Fig. Q] From a fit to the equation of state of the 
two phases we can determine the intersection point at 
fc_F^o = 25±3. This value is almost a factor two smaller 
compared to the critical density k F ro ~ 60 [2(| 28, 29[ 



using the variational method, wc optimize the width 7 of 
the stripes and their separation Ay = \y a +i—y a \- For the 
latter quantity we find that k F Ay = V^Lt is an optimal 
value corresponding to a square box having L x = L y . We 
then perform FN-DMC simulations using PBC with dif- 
ferent numbers of particles TV and we extrapolate to the 
thermodynamic limit relying onal/JV linear fit. The re- 
sults are shown in the inset of Fig. [TJ where we report the 
energy difference between the stripe and FL phase. For 
all values of k F ro in the vicinity of the freezing density 
the stripe phase is never energetically favorable compared 
to neither the FL nor the WC phase. 

In conclusion, we investigated the ground state of a 
purely repulsive system of dipolar fermions and its liquid 
to solid transition. Important extensions of this work 
concern the effect of a tilting angle, making the interac- 
tion in the 2D plane anisotropic, and the coupling to a 
second layer inducing interlayer attraction. 
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work has been supported by ERC through the QGBE 
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